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Abstract 

We propose a new functional for the charge and mass identification in AE — E 
telescopes. This functional is based on the Bethe formula, allowing safe interpolation 
or extrapolation in regions with low statistics. When applied to telescopes involving 
detectors delivering a linear response, as silicon detectors or ionization chambers, 
a good mass and charge identification is achieved. For other detectors, as caesium- 
iodide used as a final member of a telescope, a good accuracy is also obtained except 
in the low residual energy region. A good identification is however recovered if a 
non-linear energy dependence of the light output is included. 
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1 Introduction 

Stacks of detectors, called telescopes, measuring tlie energy loss and residual 
energy of cliarged particles have been used for a long time to get charge and 
mass identification, and also energy, of such particles. Several combinations 
of detectors have been used for this purpose : ionization chambers, silicon 
detectors, plastic scintillators, thallium-activated caesium-iodide scintillators 
(CsI(Tl)) read by photomultipliers or photodiodes. The identification is ob- 
tained by plotting the energy loss in one or several components of the detector 
stack versus the residual energy released in the detector in which the particle 
has stopped. In such a plot events of a given charge and mass cluster around 
identification lines. Two methods can then be implemented to extract the 
identification for each event : 
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• interactive drawing of lines on top of ridges corresponding to a given charge 
or mass, or contours around events clustered around these ridges, any event 
being identified from its distance to ridge lines or its inclusion in one contour, 

• fit of the ridge lines with a functional in which Z and A enter as parameters, 
the identification being obtained by inversion of the functional for given A£' 
and E in order to extract Z and A. 

Whereas the first method is probably more powerful and allows to face any 
situation, it suffers two main drawbacks : it doesn't deliver any extrapolation 
in regions of low statistics, and it is human-time consuming because each 
identification line has to be accurately drawn. This last aspect becomes really 
a concern when multidetectors are used, in which hundreds of such telescopes 
are involved. The second method doesn't suffer these inconveniences provided 
that the functional accurately models the data. Identification functionals have 
already been used in the past [1]. Due to their limited range of application, 
some extensions have been added later [2]. The purpose of this work is to 
propose a more complex functional based on physical grounds in order to 
allow accurate modelling of the data and safe extrapolations in regions where 
no data are present. 

The second section is devoted to the derivation of a functional based on Bethe's 
formula with a power law velocity dependence, and to its main properties. The 
third section proposes a phenomenological extension of this functional for data 
departing from the simple restriction of Bethe's formula. In the fourth section 
we will show how the introduction of a non-linear term for the amplitude of 
the light of a Csl(Tl) crystal allows a good charge identification up to Z = 50 
for telescopes involving such detectors. 



2 Basic functional 

We assume that the stopping power of a fragment of energy mass A and 
nuclear charge Z in a detecting medium takes the simple form : 



This formula can be straightly derived from Bethes's formula [4] when the 
charge state inside the stopping medium is strictly equal to the nuclear charge 

Z. Its range of validity is restricted to light ions with energy sufficient to insure 
that they are fully stripped. In particular the Bragg zone of the stopping power 
curve, where the mean charge state is no more constant, is not addressed by 
this formula. 
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We also define the integral F of / as : 

X 

F{x) = j m dt (2) 



By integration of equation (1) one obtains the range-energy relation : 

F{£/A) = (3) 

This last relation can be applied to the case of telescopes made of a first 
detector of thickness AX. The incoming fragment releases a part AE of its 
energy in this detector, and its residual energy E in the rear detector for a 
residual range X. 

When apphed to the total energy and to the residual energy, relation (3) 
reads : 

which, by elimination of X, delivers the E-AE relation : 

ae^aIf-^ 



This formula is rigorous as long as the stopping power takes the form (1). 
However to reach a practical use we need to specify the function /. We will 
choose a form close to Bethe's behaviour but simpler in order to obtain and 
analytical form of the integral inverse F"^. 



2.1 Specialisation of the function f 

In the case of Bethe's formula in the non-relativistic domain, f{£/A) is al- 
most proportional to S/A because the logarithmic term can be considered as 
constant. Therefore we will choose a power law dependence : 

f{S/A) = {81 AY (5) 
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with jj, ^ 1. This particular form of / is advantageous because it leads to an 
analytical form of (4) which reduces to : 



AE 



E"^^ + {ii + l)Z^A^' AX 



-E 



(6) 



One can notice that this relation is strictly equivalent to the functional used 
in [1] because both are based on the same hypotheses. 



2.2 Properties of the E-AE relation 



By taking the derivative of expression (6) versus E at E = one gets the 
slope at the starting point of the identification line which turns out to be 
equal to -1 independently of Z and A. For signals delivered by real detectors 
this slope includes the ratio of electronic gains. However, in case of detectors 
exhibiting a linear response as sihcon detectors, this constancy of the starting 
slope is generally achieved. Conversely the verification of this property is a 
good indication of the validity of approximations brought by (1). 

By setting = in (6) one obtains the series of crossing points of identification 
lines with the energy loss axis : 

AEo = [{^i+l)AX]^Zl^Al^ (7) 



In particular for = 1, close to reality, these starting energy losses are pro- 
portional to Zy/A. 

At high energy the second term between [. . .] in equation (6) becomes small, 
compared to the first one, and an expansion to first order can be performed : 

which is proportional to the stopping power, as expected. 



2.3 Power law functional 



If we assume that the detector response is linear, and that data comply to 
relation (6), then they should be fitted by the following function : 

1 

- gE (9) 
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where the parameters are g (ratio of electronic gains), A which includes the 
thickness of the first detector and /i which is close to 1. Furthermore E and AE 
are net signals from which pedestals have been subtracted. We then obtain a 
3-parameter formula, or a 5-parameter formula if we also fit on the AE and 
E pedestals. 

To guarantee a convergence of the fit in all situations one should provide rea- 
sonable starting values and impose constraints on the range of the parameters. 
In particular fj, should lie between 0.5 and 1.5 with a starting value /Iq — 1, the 
initial value Aq of A should be determined from the starting point of an identi- 
fication line, and go should be obtained from the slope at the same point. The 
reader is referred to the appendix for technicalities related to the fit procedure. 

Once parameters have been determined by the fitting procedure, the func- 
tion (9) can be used to extract Z and eventually A associated to any AE — E 
pair. This is achieved by an analytical inversion of equation (9) which deliv- 
ers the quantity Z'^A'^ as early observed [1]. If one is only interested in the 
charge identification A is set dependent on Z, the simplest prescription being 
A — 2Z, and equation (9) is solved for Z and the solution is projected onto the 
nearest integer number Z^. Furthermore if the mass identification is required, 
the value Zi is taken for Z and the equation is solved for A. 

As an example figure 1 shows the AE — E map obtained with a calibration 
module of INDRA, made of a silicon detector of thickness 70 /xm and a lithium- 
diffused silicon detector 2 mm thick, placed in front of a CsI(Tl) crystal used 
as a veto removing all particles punched through the silicon-lithium detector. 
The 8 dashed lines are the reference lines, from hydrogen to beryllium, used 
as references for the 5-paramctcr fit. The result of the fit is shown on figure 2 
which displays the particle identification PI = Zi + 0.2{A — 2Zi) versus the en- 
ergy released in the silicon-lithium detector. In can be verified that distorsions 
are very small and that the extrapolation toward boron an carbon elements 
is correct, indicating that the power law formula models correctly the data in 
this range. 

However it has already been noticed in the past [2,3] that for a wider Z range, 
the above formula must be extended. The stopping power may depart from 
the behaviour ruled by relations (1) and (5), moreover the response of some 
detectors is not completely linear : pulse height defect in silicon detectors, light 
response of Csl(Tl) detectors. . . Figure 3 shows the AE — E map obtained 
with a telescope made of 2 silicon detectors, of thicknesses equal to 150 iim 
and 750 /im, in the case of a ^°^Pb beam at 29 MeV/n [8]. Thick dashed lines 
drawn on top of ridges associated to Z = 8,13,17,22,28,37 are used as a 
reference for the 5-parameter fit with relation (9) which generates the thin 
lines. It turns out that the functional hardly reproduces the data in the whole 
range : not only its Z dependence is approximative, but also the shape of the 
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Fig. 1. AE — E map from a silicon-silicon telescope. The 8 dashed lines, from 
hydrogen to beryllium, are the reference lines used for the 5-parameter fit of equa- 
tion (9), which in turn generates the thin continuous lines. For Z = 6, 7, 8, A = 2Z 
is assumed. 

energy dependence for a given line. 

Therefore an enrichment of (9) is highly desired in order to benefit from more 
degrees of freedom. The extension proposed hereafter is not based on theoret- 
ical considerations but rather on phenomenological requirements. 



3 Extended functional 



One can see from equation (9) that /i rules the ordinates of starting points of 
identification lines. Once ^ has been fitted to reproduce these data, the rapidity 

2 M 

of the transition from the low energy domain dominated hj \Z Ai^+^ to the 
high energy region, proportional to the stopping power, is fixed. This could 
lead to a too strong constraint when applied to real data. One could introduce 
an additional parameter for the Z exponent but this would still link the Z 
dependencies in the low and high energy regions. As regards the A exponent, 
it can't be touched because it insures the right asymptotic behaviour at high 
energy, toward the stopping power. 
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Fig. 2. Plot of the particle identification variable defined as PI = Zi + 0.2(A — 2Zj) 
where Zi is the nearest integer deduced from the functional (9) setting A = 2Z, and 
A is the mass which solves (9) when Z = Z^. 

The ideal solution should meet the following prescriptions : independent vari- 
ations in Z and A in the low and high energy regimes, possibility of tun- 
ing separately the rapidity of the transition between these regimes, correct 
asymptotic behaviour at high energy. This is achieved by the addition of new 
phenomenomogical term. 



We define the extended functional as : 

.fJ.+ U+l 



{gEY+^+^ + (XZ'^A'^Y + iZ^A^'{gEy 



gE (10) 



which reduces to the basic functional (9) when one substitutes iJ, for 11 + u and 
sets a — 2/{ijl+1), P — + and ^ = 0, or even more simply if one sets : 
u — and A = 0. 

The g parameter still represents the slope of identification hues at their start- 
ing points (at E — 0), equal for all lines as for the power law formula. The 
energy loss value at the starting points writes : AEq = XZ^'A'^, showing that 
this series of points essentially determines the A, a and /? parameters. 

At higher energies an expansion of (10) to first order can be carried out. 
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Fig. 3. AE — E map obtained with a 150 ;um-750 jim. silicon telescope for the 
208p|-,_)_natg- j-eaction at 29 MeV/n [8]. More populated areas are indicated by clearer 
zones. Thick dashed hnes drawn for Z = 8, 13, 17, 22, 28, 37 are used as a reference 
for the 5-parameter fit with relation (9) which generates the thin lines. 

leading to : 



X 



+ {gE/Ay 



X 



_^ V ) 



{gE/AY 



One can notice that at high energy the dE/ dX is recovered if ji stays in 
the vicinity of 1, and that ^ governs the amount of energy loss in this region. 
Moreover one benefits from the second term inside the square brackets which 
allows, through the v parameter, to change the shape of curves in the region 
of intermediate energy. 

The extended formula is a 7-parameter functional (A, a, yU, z/, ^ and g) or 
a 9-parameter one if both pedestals are also considered as parameters of the 
fit. The convergency of the fit is always insured only if good starting values 
are supplied and if constraints are applied to the range of parameters. More 
precisely one should adopt : aQ — — 1, I3q — 0.5, vq — 1, Aq determined 
from the starting point of an identification line, g^ from its starting slope, 
and ^0 from one point taken at high energy. The ranges should be restricted 
to : 0.5 < a < 1.5, 0.2 < /3 < 1, 0.2 < < 1.5, 0.1 < z/ < 4, 
Ao/4 < A < 4Ao, 5'o/4 < g < 4:go, < ^. In the case where no mass 
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Fig. 4. Same as figure 3 but thin lines result now from the fit based on the 8-pa- 
rameter extended functional (10), assuming A = 2Z. 

identification is given, the number of parameters is decremented by one : j3 is 
kept equal to 0.5 and A is taken as a function of Z, the simplest relationship 
being A = 2Z. 

Unlike the power law formula, the extended formula can't be analytically 
inverted to extract Z and A from a AE — E pair. This can only be done nu- 
merically and the fastest algorithm is probably the Newton-Raphson method 
which however needs the derivatives of (10) with respect to Z and A. 

The quality of a fit based on this extended formula with 8 parameters, as- 
suming A = 2Z and (3 arbitrarily set, is shown on figure 4 where it appears 
that the shape of the curves is nicely reproduced, and also the interpolation 
in regions where no reference line guides the fit procedure. Even the extrapo- 
lation toward higher Z's is satisfactory, having in mind that in this region the 
simple hypothesis A = 2Z is unlikely to hold. Figure 5 shows that only small 
distorsions remain when the identified Z for each event is plotted versus the 
residual energy. 

We can also check the ability of this functional to reproduce the data for 
telescopes using a Csl(Tl) crystal as a stopping detector, for which the light 
response is not linear with the deposited energy and depends on the nature of 
the particle. For this purpose we worked on the data collected by a module of 
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Fig. 5. For the same data as in figures 3 and 4, the Z extracted from relation (10) 
is plotted versus the residual energy to exhibit possible distorsions. 




Light Light 

Fig. 6. AE — h plot for a silicon-CsI(Tl) telescope, the abscissa is the total light 

output reconstructed from the fast and slow components of the CsI(Tl). The thick 
dashed lines follow the ridges of charges Z = 4, 6, 12, 22, 30, 40 and the thin lines 
result from the 8-parameter fit based on the extended functional (10) with A = 2Z. 
The right part is a zoom applied to the low light region of the left map, including the 
line Z = 30. The non-linearity of the light response shows up by the high curvature 
of data, not predicted by the functional. 
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the second ring of the muhidetector INDRA, made of a 300 fim sihcon planar 
detector followed by a CsI(Tl) of depth equal to 14 cm [9]. As the CsI(Tl) 
signal was only partially integrated on 2 time intervals for the delivering of 
the so-called slow and fast components, we didn't have at disposal a true 
measurement of the total light signal, integrated over the duration of the 
pulse. We reconstructed this total light signal h by a combination of slow 
and fast components following a procedure described in [7]. Figure 6 is an 
illustration of a AE — h plot collected with a tin beam. The thick dashed lines 
are drawn on top of charges Z = 4, 6, 12, 22, 30, 40 and the thin lines result 
from the 8-parameter fit with A — 2Z. Despite the non-linearity of the light 
response and its dependence with Z the fit is in good agreement with the data, 
showing that the parameters of the functional are able to partially compensate 
for these effects. However by looking carefully at the low residual energy region, 
as displayed in the right part of figure 6, one can see that the curvature of 
data is much more pronounced than predicted by the functional. This effect is 
due to the highly non-linear light response which is almost quadratic in this 
region. Any improvemevent of the accuracy of generated identification lines in 
this region needs an explicit inclusion of such a non-linearity. 



4 Light response of caesium-iodide crystal 

The idea for a new extension of functionals (9) and (10) is to substitute for 
E its evaluation from the light output, in order to derive direct dependencies 
of the energy loss with the light output, the most part of non-linearities being 
carried by the relationship between light and energy. 

We will take for the light response of CsI(Tl) crystals a dependence deduced 
from the Birks formula [5]. If one denotes h the net light output after sub- 
traction of the pedestal, its dependence with the energy deposition takes the 
form [6] : 



and ?7 is a new parameter. No multiplicative constant is needed because the 
scaling factor is already included in the g parameter of the functional, this 
is equivalent to consider the energy and the light output as measured with 
the same unit. This formula is not the most accurate [7] and particularly it 
discards the loss of quenching due to the delta electrons. We take it merely 
as a reasonable way of accounting for the non-linearities at low energy and 
high Z, hoping that the other parameters of the functional will be able to 
compensate its partial inadequacy. 




where 




(11) 



11 



However one difficulty arises from the fact that equation (11) cannot be ana- 
lytically inverted to express £^ as a function of as required for its insertion 
into the functionals (9) or (10). This is the reason why we adopted a slightly 
different formula for the inverse relation, imposing the same asymptotic be- 
haviours at low and high energy. Prom (11) we derive for the low and high 
energy regimes : 



^2 

E — ^ ^ h — > — 
E — ^oo =^ h — >E-p\n(^^ 

which exhibits the quadratic dependence at low energy. By inverting these 
relations one gets : 



h 



E 



2ph 



h 



oo 



E — > h + pln 



Now we seek for a formula which complies with the above limiting behaviours, 
and we can check; that : 



E 



\ 



h? + 2ph 



1 + In 1 + 



h 



[12) 



fulfills this condition. Obviously relation (12) is not the strict inverse of (11), 
but it departs from it by 7.4 % only at maximum, and both become identical 
at low and high energy. 

If we use this relation in combination with (10), we obtain now a functional 
with 8 parameters (A, a, /?, /x, z/, ^, g and rj), and 2 eventually additional 
parameters if pedestals in h and A£J are also adjusted. 

Figure 7 displays the calculated lines resulting from a 9-parameter fit in which 

/3 has been kept constant as the charge identification only is given, and A = 2Z 
is assumed. It can be seen that a satisfactory agreement is obtained for the 
shape of the curves. Particularly one can notice that the high curvature of 
the lines in the vicinity of the energy loss axis is well reproduced by the 
introduction of the light response and its quadratic dependence at low energy. 
Furthermore a stable extrapolation is deduced for the charges higher than 40, 
the highest line used as a reference for the fit. 

The quality of the charge identification can be checked on figure 8 where the 
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Fig. 7. Same as figure 6 but the thin hues result now from a 9-parametcr fit 
combining relations (10) and (12). The curvature of data at low energy is now 
reproduced. 
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Light 

Fig. 8. For the same data as in figures 6 and 7, the Z extracted from relations (10) 
and (12) is plotted versus the residual energy to exhibit possible distorsions. 

distorsions of the identified Z versus the hght output appear to be small even 
for a wide range of charges. The apparent loss of resolution at low Z's simply 
reflects the coding granularity, the channel size in energy loss becoming too 
large for these low charges. 
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5 Conclusion 



Starting from the usual power law formula we propose an extended functional 
for the particle identification in AE" — E telescopes, allowing a good mass 
and charge identification with silicon detectors on a wide range in energy 
and charge. Its domain of applicability remains however limited to situations 
where the Bragg region plays a minor role, in particular this excludes high 
Z's detected in thin AE detectors. This functional can also be applied to spe- 
cific cases where the E signal is no longer linear with the energy deposition, 
as for example for silicon-Csl(Tl) telescopes in which the fast component of 
the Csl(Tl) is taken as a measure of the residual energy. This procedure has 
been applied to the calibrating modules in the INDRA multidetector, how- 
ever it appears to be accurate only on a range of a few charges. By using the 
total hght output of CsI(Tl), either by integrating the whole light signal or 
by estimating it from the partial integrations associated to the slow and fast 
components, the agreement between data and the extended formula is highly 
improved on a large range of Z. Moreover the discrepancies which remain at 
low residual energy are removed by the explicit introduction of the non-linear 
light response of the CsI(Tl) crystal, and a good quality of the identification 
is obtained on a wide range in charge and energy. This new functional opens 
the possibility of identifying particles in silicon or silicon-Csl(Tl) telescopes by 
defining only a limited number of reference points, and allows to rely on the 
deduced interpolation and even on the extrapolation in low statistics regions. 
It should make easier the usual identification task of particle identification in 
the case of charged particle multidetectors. The author expresses his aknowl- 
edgement to M. Morjean and to the INDRA collaboration, for having provided 
their data before publication. 



A Appendix 

The fit algorithm relies on the local quadratic expansion of the x^- At each 
step a calculation is made for x^, its gradient, and an estimation of the hessian 
matrix obtained by neglecting the second derivatives of the functional. This 
approximation is powerful because it always delivers a definite positive matrix 
which, in the vicinity of the minimum, becomes very close to the true hessian 
matrix. With such a prescription it is necessary to compute only the func- 
tional value and its gradient at each data point. The minimization proceeds 
iteratively and at each step the status of constraints is checked. 

In order to insure a correct convergency 2 conditions have to be met : 

• good starting values of the parameters have to be supplied as already ex- 
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plained, 

• for the minimization procedure the parameters have to be scaled in order 
to be of the order of unity. 

For a practical use, the formulas giving the derivatives of the functional (10) 
against parameters are written hereafter. By denoting : 



G 



the derivatives can be expressed as : 



dAE 
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da 
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dp 
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G\n{G) \ 
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The derivative dAE/dE is involved in the E'-pedestal dependence. It is also 
useful when E is expressed from the hght output using (12). The dAE/dZ 
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and dAE/dA derivatives are ingredients of the Newton- Raphson method for 
the determination of Z and A for a given AE — E pair whence the parameters 
have been determined. 

Regarding the hght response (12), the useful derivatives are : 

dE__l^{ 1 l + ln(l + /i/p) | 

'dh~E\ ^ 1 + hfp J 

dE h \ 1 , ,^ ^ , ,1 
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